# Script for "Direct Experience with Poor Working Conditions and Partisan Labor Policy Preferences" 

# OVERVIEW
# Steps:
  # 1 Download data file and add to working directory
  # 2 Run below script to read in and clean data, run models, and produce tables/figures
  # 3 View outputs as R objects: table1, table2, fig1, fig2, fig3, fig4


# Packages ----
library(haven)
library(dplyr)
library(tidyr)
library(readr)
library(ggplot2)
library(stargazer)


# Functions -----

re_workplace <- function(x) {
  x <- ifelse(x == 1, "Big-box superstore",
              ifelse(x == 2, "Department store",
                     ifelse(x == 3, "Retail store",
                            ifelse(x == 4, "Grocery store",
                                   ifelse(x == 5, "Restaurant",
                                          ifelse(x == 6, "Fast food",
                                                 ifelse(x == 7, "Coffee shop",
                                                        ifelse(x == 8, "Hotel or motel",
                                                               ifelse(x == 9, "Warehouse",
                                                                      ifelse(x == 10, "Fulfillment center",
                                                                             ifelse(x == 11, "Delivery vehicle",
                                                                                    ifelse(x == 12, "Convenience store",
                                                                                           ifelse(x == 13, "Drug store",
                                                                                                  ifelse(x == 14, "Corporate office",
                                                                                                         ifelse(x == 15, "Call center",
                                                                                                                ifelse(x == 16, "Other", x))))))))))))))))
  x
}

# Data ------

s <- read_csv("poq_data.csv")

# Clean -----

sc <- s %>%
  mutate(
    pol_unionize = as.numeric(unionpolicy == 1),
    s_less_twoweeks = ifelse(notice %in% c(1, 2), 1, 0),
    s_keep_open = as.numeric(keepschedopen == 1),
    s_on_call = as.numeric(oncall == 1),
    s_cancel_shift = as.numeric(cancelshift == 1),
    s_time_change = as.numeric(timing == 1),
    s_clop = as.numeric(clopening == 1),
    dem_pid3 = ifelse(party_id == 1 | (party_id %in% 3:4 & party_id_lean == 1), "Democrat",
                      ifelse(party_id == 2 | (party_id %in% 3:4 & party_id_lean == 2), "Republican", "Independent")),
    dem_educ = as.numeric(educ),
    dem_income = as.numeric(hhincome),
    dem_age = ifelse(age_clean > 7, median(s$age_clean), age_clean),
    dem_union = as.numeric(union),
    dem_ideo = as.numeric(ideology),
    dem_married = as.numeric(cohabstatus == 1),
    dem_ownhome = as.numeric(ownrenthome == 1),
    dem_yrs_res = as.numeric(residence),
    dem_female = as.numeric(gender == 2),
    dem_black = as.numeric(black == 1),
    dem_hisp = as.numeric(hisp == 1),
    dem_asian = as.numeric(aian == 1),
    dem_other = as.numeric(oth == 1),
    workplace = as.numeric(workplace),
    workplace_name= re_workplace(workplace),
    work_nosick =  as.numeric(benefits_paidsick == 0),
    work_novacation = as.numeric(benefits_paidvacation == 0),
    work_noleave = as.numeric(benefits_paidleave == 0),
    work_nohealth = as.numeric(benefits_health == 0),
    work_nodental = as.numeric(benefits_dental == 0),
    work_badscale = (s_less_twoweeks + s_keep_open + s_on_call + s_cancel_shift + s_time_change + s_clop +
                       work_nosick + work_novacation + work_noleave + work_nohealth + work_nodental)/11,
    work_employer = as.numeric(q1_employer),
    work_state = as.numeric(statelist_harm)) %>%
  dplyr::select(starts_with("s_"),
                starts_with("pol_"),
                starts_with("dem_"),
                starts_with("work_"),
                workplace, workplace_name)


sc <- select(sc, pol_unionize, work_badscale, dem_pid3, dem_ideo, dem_educ,
             dem_income, dem_union, dem_age, dem_female,
             dem_married, dem_yrs_res, dem_ownhome,
             dem_black, dem_hisp, dem_asian, dem_other,
             workplace_name, work_state, work_employer,
             s_less_twoweeks, s_keep_open, s_on_call, s_cancel_shift, s_time_change, s_clop,
             work_nosick, work_novacation, work_noleave, work_nohealth, work_nodental)
sc <- sc[complete.cases(sc),]

# Figure 1: Job attributes among Democratic and Republican workers ------
bji <- select(sc, dem_pid3, s_less_twoweeks, s_keep_open, s_on_call, s_cancel_shift, s_time_change, s_clop,
              work_nosick, work_novacation, work_noleave, work_nohealth, work_nodental)

fig1 <- bji %>%
  filter(dem_pid3 != "Independent") %>%
  group_by(dem_pid3) %>%
  summarise_all(mean, na.rm=T) %>%
  pivot_longer(cols = 2:12) %>%
  ggplot(aes(x = reorder(name, value),
             y = value,
             fill = dem_pid3)) +
  geom_col(position = "dodge") +
  coord_flip() +
  scale_fill_grey() +
  scale_y_continuous(breaks = c(0, .10, .20, .30, .40, .50, .60, .70),
                     labels = c("0%", "10%", "20%", "30%", "40%", "50%", "60%", "70%")) +
  scale_x_discrete(labels = c("Had canceled shift",
                              c("Had to be on-call within hours of work",
                                "No health insurance",
                                "No vacation",
                                "No dental",
                                "Had to work 'clopening'",
                                "No sick leave",
                                "Received schedule less than two weeks in advance",
                                "Had timing or length of shift changed",
                                "No paid parental leave",
                                "Had to keep schedule open to be available for job"))) +
  theme_minimal() +
  labs(y = "",
       x = "",
       caption="Note: Percentages show % of workers, by party,\nwho indicated the attribute was present in their job") +
  theme(legend.position = "bottom",
        legend.title = element_blank(),
        text = element_text(family = "Times"))

# Figure 2: Distribution of PWCI, full sample -------
fig2 <- sc %>%
  filter(!is.na(dem_pid3)) %>%
  ggplot(aes(x = work_badscale)) +
  geom_histogram(bins = 12) +
  theme_minimal() +
  labs(x = "Poor Working Conditions Index (PWCI)",
       y = "n",
       caption = "bins=12") +
  theme(text = element_text(family = "Times"))

# Figure 3: PWCI density among Democratic and Republican workers -------
fig3 <- sc %>%
  filter(dem_pid3 != "Independent") %>%
  ggplot(aes(x = work_badscale,
             color = dem_pid3)) +
  geom_density(alpha = 0.1) +
  theme_minimal() +
  labs(x = "Poor Working Conditions Index (PWCI)") +
  scale_color_grey() +
  theme(legend.position = "bottom",
        legend.title = element_blank())

# Table 1: models -------
a1 <- glm(pol_unionize ~ work_badscale, data = sc, family = binomial("logit"))
a2 <- glm(pol_unionize ~ work_badscale + dem_pid3 + dem_ideo, data = sc, family = binomial("logit"))
a3 <- glm(pol_unionize ~ work_badscale + dem_pid3 + dem_ideo + dem_educ + dem_income + dem_union + dem_age + dem_female +
            dem_married + dem_yrs_res + dem_ownhome + dem_black + dem_hisp + dem_asian + dem_other, data = sc, family = binomial("logit"))
a4 <- glm(pol_unionize ~ work_badscale + dem_pid3 + dem_ideo + dem_educ + dem_income + dem_union + dem_age + dem_female +
            dem_married + dem_yrs_res + dem_ownhome + dem_black + dem_hisp + dem_asian + dem_other + factor(work_state), data = sc, family = binomial("logit"))

a5 <- glm(pol_unionize ~ work_badscale + dem_pid3 + dem_ideo + dem_educ + dem_income + dem_union + dem_age + dem_female +
            dem_married + dem_yrs_res + dem_ownhome + dem_black + dem_hisp + dem_asian + dem_other + workplace_name, data = sc, family = binomial("logit"))

a6 <- glm(pol_unionize ~ work_badscale + dem_pid3 + dem_ideo + dem_educ + dem_income + dem_union + dem_age + dem_female +
            dem_married + dem_yrs_res + dem_ownhome + dem_black + dem_hisp + dem_asian + dem_other + factor(work_employer), data = sc, family = binomial("logit"))
a7 <- glm(pol_unionize ~ work_badscale + dem_pid3 + dem_ideo + dem_educ + dem_income + dem_union + dem_age + dem_female +
            dem_married + dem_yrs_res + dem_ownhome + dem_black + dem_hisp + dem_asian + dem_other + workplace_name + factor(work_state) + factor(work_employer), data = sc, family = binomial("logit"))

# Table 1: table -----
varlabels1 <- c("Poor Working Cond. Index","Independent","Republican",
                "Ideology","Education","Income","Union member","Age",
                "Female","Married","Residential permanence","Own home",
                "Black","Latino","Asian","Other race","Intercept")

table1 <- stargazer(list(a1, a2, a3, a4, a5, a6, a7),
          type = "text",
          header = F,
          font.size = "scriptsize",
          omit = "factor|workplace_",
          add.lines = list(c("State Fixed Effects", "No", "No", "No", "Yes", "No", "No", "Yes"),
                           c("Workplace Type Fixed Effects", "No", "No", "No", "No", "Yes", "No", "Yes"),
                           c("Employer Fixed Effects", "No", "No", "No", "No", "No", "Yes", "Yes")),
          covariate.labels = varlabels1,
          digits = 3,
          report = "vcp",
          omit.table.layout="n")



# Table 2: models -------
aa1 <- glm(pol_unionize ~ work_badscale, data = sc, family = binomial("logit"))
aa2 <- glm(pol_unionize ~ work_badscale*dem_pid3 + dem_ideo, data = sc, family = binomial("logit"))
aa3 <- glm(pol_unionize ~ work_badscale*dem_pid3 + dem_ideo + dem_educ + dem_income + dem_union + dem_age + dem_female +
             dem_married + dem_yrs_res + dem_ownhome + dem_black + dem_hisp + dem_asian + dem_other, data = sc, family = binomial("logit"))
aa4 <- glm(pol_unionize ~ work_badscale*dem_pid3 + dem_ideo + dem_educ + dem_income + dem_union + dem_age + dem_female +
             dem_married + dem_yrs_res + dem_ownhome + dem_black + dem_hisp + dem_asian + dem_other + factor(work_state), data = sc, family = binomial("logit"))

aa5 <- glm(pol_unionize ~ work_badscale*dem_pid3 + dem_ideo + dem_educ + dem_income + dem_union + dem_age + dem_female +
             dem_married + dem_yrs_res + dem_ownhome + dem_black + dem_hisp + dem_asian + dem_other + workplace_name, data = sc, family = binomial("logit"))

aa6 <- glm(pol_unionize ~ work_badscale*dem_pid3 + dem_ideo + dem_educ + dem_income + dem_union + dem_age + dem_female +
             dem_married + dem_yrs_res + dem_ownhome + dem_black + dem_hisp + dem_asian + dem_other + factor(work_employer), data = sc, family = binomial("logit"))
aa7 <- glm(pol_unionize ~ work_badscale*dem_pid3 + dem_ideo + dem_educ + dem_income + dem_union + dem_age + dem_female +
             dem_married + dem_yrs_res + dem_ownhome + dem_black + dem_hisp + dem_asian + dem_other + workplace_name + factor(work_state) + factor(work_employer), data = sc, family = binomial("logit"))

# Table 2: table ------
varlabels2 <- c("PWCI","PWCI x Independent", "PWCI x Republican",
                "Independent","Republican",
                "Ideology","Education","Income","Union member","Age",
                "Female","Married","Residential permanence","Own home",
                "Black","Latino","Asian","Other race","Intercept")


table2 <- stargazer(list(aa1, aa2, aa3, aa4, aa5, aa6, aa7),
          type = "text",
          header = F,
          font.size = "scriptsize",
          omit = "factor|workplace_",
          order = c("work_badscale","dem_pid3Independent","dem_pid3Republican",
                    "work_badscale:dem_pid3Independent","work_badscale:dem_pid3Republican",
                    "dem_ideo",
                    "dem_educ",
                    "dem_income",
                    "dem_union",
                    "dem_age",
                    "dem_female",
                    "dem_married",
                    "dem_yrs_res",
                    "dem_ownhome",
                    "dem_black",
                    "dem_hisp",
                    "dem_asian",
                    "dem_other",
                    "Constant"),
          covariate.labels = varlabels2,
          add.lines = list(c("State Fixed Effects", "No", "No", "No", "Yes", "No", "No", "Yes"),
                           c("Workplace Type Fixed Effects", "No", "No", "No", "No", "Yes", "No", "Yes"),
                           c("Employer Fixed Effects", "No", "No", "No", "No", "No", "Yes", "Yes")),
          digits = 3,
          report = "vcp",
          omit.table.layout="n")



# Figure 4: predicted probability plot -----
ppaa1 <- ggeffects::ggpredict(aa3, c("work_badscale", "dem_pid3"), typical = "median")

fig4 <- ppaa1 %>%
  filter(group != "Independent") %>%
  ggplot(aes(x = x,
             y = predicted,
             linetype = group)) +
  geom_line() +
  geom_ribbon(aes(ymax = conf.high,
                  ymin = conf.low),
              alpha = 0.2) +
  theme_minimal() +
  labs(x = "Poor Working Conditions Index (PWCI)",
       y = "Predicted probability") +
  theme(legend.title = element_blank(),
        text = element_text(family = "Times"))




